##  Balance test for studies 1-3 for sponsorship paper
rm(list = ls())


setwd("~/Dropbox/Current projects/kern-pietryka-crabtree/data/")
load("Experiment 2 intermediate results.RData")


##  Balance test for Trump vote share in 2016
##  We find fips code and add vote share data
votes <- read.csv("US_County_Level_Presidential_Results_12-16.csv", header = TRUE)

##  this version of the code uses a different file from 
##  https://github.com/tonmcg/US_County_Level_Election_Results_08-16
##  than the previous version
##  results are of course identical


##  reformat fips codes
sel <- nchar(votes$combined_fips) == 4
table(sel)
votes$combined_fips[sel] <- paste("0", votes$combined_fips[sel], sep = "")
table(nchar(votes$combined_fips))


dat$Trump <- NA

m <- 1
for(m in 1:nrow(dat))	{

	temp <- latlong2fips(latitude = as.numeric(as.character(dat$LocationLatitude[m])),
						 longitude = as.numeric(as.character(dat$LocationLongitude[m])))

if(is.null(temp))	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }
if(is.na(temp)) 	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }
if(temp == "NULL")	{	cat(m,as.numeric(as.character(dat$LocationLatitude[m])),
						as.numeric(as.character(dat$LocationLongitude[m])),"\n"); next() }

	sel <- which(votes$combined_fips == temp)
	stopifnot(length(sel) == 1)

	dat$Trump[m] <- votes$per_gop_2016[sel]

	if(m %% 100 == 0)	{ cat(m, "\n") }

						}

##  consistency check
sum(is.na(dat$Trump))
range(dat$Trump)


##  RI for balance check using difference in means from OLS model
set.seed(1)
t.bal.rand <- matrix(NA,1000000,3)

m <- 1
for(m in 1:1000000)	{

	dat$Tr.r <- sample(c(1,2,3), nrow(dat), replace = TRUE)
	t.bal.rand[m,] <- lm(dat$Trump ~ as.factor(dat$Tr.r))$coef

	if(m %% 1000 == 0)	{ cat(m, "\n") }

					}

Omega.bal <- solve(var(t.bal.rand[,2:3]))

t.bal <- lm(dat$Trump ~ as.factor(dat$Tr))$coef
t.bal <- t(t.bal[2:3]^2) %*% Omega.bal %*% t(t(t.bal[2:3]^2))


out.bal <- rep(NA,1000000)
m <- 1
for(m in 1:nrow(t.bal.rand))	{

	out.bal[m] <- t.bal.rand[m,2:3]^2 %*% Omega.bal %*% t(t(t.bal.rand[m,2:3]^2))

								}

p.bal <- round(mean(as.vector(t.bal) < out.bal),2)
p.bal















.


